Неделя один

Обозначения

  • Одна зависимая (объясняемая) переменная: y

  • Несколько регрессоров (предикторов, объясняющих переменных): x, z…

  • По каждой переменной n наблюдений: \(y_1, y_2...y_n\)

Модель — это формуля для объясняемой переменной.

Пример

Возьмём например данные по машинам 1920-х годов. Тут видимая линейная взаимосвязь.

Данные по машинам 1920-х гдов

Данные по машинам 1920-х гдов

Модель для этих данных может иметь вид \(y_i = \beta_1 + \beta_2x_i + \varepsilon_i\)

Что здесь есть что? У нас есть наблюдаемые переменные,

  • y — это длина тормозного пути

  • x — это скорость, с которой ехала машина.

  • Есть неизвестные коэффициенты $ _1, _2$

  • случайная составляющая, ошибка

То есть, \(\beta_2\) показывает — насколько увеличивается тормозной путь, если машина разгонится на один лишний километр в час. И есть некая случайная составляющая \(ε_i\), это может быть все что угодно:

  • водитель по-другому нажимал на тормоз,
  • что-то там на дороге было другое,

то есть это та часть, по которой у нас нет возможности предсказать, но, тем не менее, вот эта случайная ошибка она входит в y.

План действий

  • Придумать адекватную модель

  • Получить оценки неизвестных параметров \(\widehat\beta_1, \widehat\beta_2\)

  • Спрогнозировать, заменив неизвестные параметры на оценки \(\widehat y_i = \widehat\beta_1 + \widehat\beta_2 x_i\)

Метод наименьших квадратов

Как найти \(\beta_1, \beta_2\)? Собственно методом наименьших квадратов.

Если мы придумали какие-то оценки, \(\beta_1, \beta_2\) то соответственно возникает такое понятие как ошибка прогноза

\[\widehat \varepsilon_i = y_i - \widehat y_i\]

Есть суммарная ошибка, чтобы суммарная ошибка не занулялась одна в плюс, другая в минус, не компенсировали друг друга, мы возведем в квадрат. И посчитаем сумму квадратов ошибок прогноза, то есть сумма \(\widehat \varepsilon_i ^2\)

\[Q(\widehat\beta_1, \widehat\beta_2) = \sum_{i=1}^{n} \widehat\varepsilon^{2}_{i} = \sum_{i=1}^{n} (y_i - \widehat y_i)^2 \]

Суть метода МНК

Возьмите в качестве оценок такие коэффициенты \(\widehat\beta_1, \widehat\beta_2\) при которых сумма квадратов ошибок прогноза будет минимальна

Случай для одного регрессора

Всё сводится к тому, что в модели \[M_1: y_i = \beta + \varepsilon_i\]

\[\widehat\beta = \bar y\]

Случай для одного регрессора

Случай для одного регрессора

Случай для двух регрессоров

Тут несколько сложнее, подробности на скриншоте

Случай для двух регрессоров

Случай для двух регрессоров

Случай для множества регрессоров

Что мы получили:

Случай для двух регрессоров

Случай для двух регрессоров

Ещё раз пробежимся по терминам

Случай для двух регрессоров продолжение

Случай для двух регрессоров продолжение

Графическое представление

Изобразим на графике, всё то, о чем мы только что проговорили.

Графический вывод

Графический вывод

В частности \(\widehat\beta_1\), т.е. точка на оси ординат от пересечении с прямой регрессии, называют пересечением или Intercept-ом

Метод наименьших квадратов подбирает прямую так, чтобы суммарные расстояния от точек до прямой были минимальными.

Графический вывод продолжение

Графический вывод продолжение

Множественные регрессоры

Случай множественных регрессоров принципиально не отличается от двух регрессоров. Поэтому рассмотрим на примере трёх предикторов.

Множественные регрессоры

Множественные регрессоры

Суммы квадратов

Суммы квадратов

Суммы квадратов

Что есть что:

  • RSS (Residual Sum of Squares) — этот показатель меряет, насколько велики эпсилон с крышкой, насколько они далеко лежат от нуля.

  • TSS (Total Sum of Squares) — этот показатель меряет, насколько каждый из \(y_i\) не похож на \(y\) среднее. Если \(y_i\) далеко лежит от \(y\) среднее, соответственно, это слагаемое в сумме квадратов будет большим. И вся сумма квадратов будет большой.

  • ESS (Explained Sum of Squares) — она показывает, насколько прогнозное значение \(y_i\) с крышкой далеко легло от \(y\) среднего.

Лекбез по линейной алгебре

Язык эконометрики во многом — это язык линейной алгебры, нужно его знать.

Обозначения

Обозначения

Обозначения

  • Маленькой буквой \(y\) мы будем обозначать вектор, то есть столбик из всех игреков, записанных друг под другом: у_1, у_2 и так далее, у_n.

  • Ну, соответственно, \(х\) маленькое — это все иксы: х_1, х_2 и так далее, х_n, записанные друг под другом.

  • Аналогично \(\widehat\beta\).

  • И еще введем обозначение: единичку со стрелочкой. Это будет вектор-столбец, то есть столбик из единичек в количестве n штук, потому что у нас в модели будет n наблюдений.

Тогда для нашей модели

\[\widehat y = \widehat\beta_1 \cdot \overrightarrow 1 + \widehat\beta_2 \cdot x + \widehat\beta_3 \cdot z\]

  • Большая буква X — это все регрессоры, помещенные в одну большую табличку чисел, которые называются матрицей.
Таблица регрессоров

Таблица регрессоров

  • Рассмотрим ещё такое понятие как длина вектора
Длина вектора

Длина вектора

Где у нас бывают длины векторов:

Пример длины вектора

Пример длины вектора

Есть одно замечательное применение. Если скалярное произведение векторов равно нулю, значит они перпендикулярны.

Скалярное произведение векторов

Скалярное произведение векторов

x <- c(1,2, -3)
y <- c(-6,0,-2)

sum(x * y)
## [1] 0

Линейная модель регрессии в n-мерном пространстве

n-мерное простарнство

n-мерное простарнство

Есть вектор y, есть вектор из одних 1, я продолжаю этот вектор до прямой и оказывается, что

  • в простой модели без регрессоров спроецировав вектор y на эту прямую, я получу вектор y с крышкой — предсказанные значения.

Соответственно, мы получили попутно еще один факт: если любой вектор проецировать на вектор из 1, то получится вектор средних значений

Геометрия множественной регрессии

Напомню, что мы вывели условие первого порядка

Геометрическая интерпретация условий первого порядка

Геометрическая интерпретация условий первого порядка

Это означает, что вектора перпендикулярны.

Облачко это все те вектора, которые можно получить, складывая с некоторыми весами вектор x, вектор z и вектор из единичек. Это гиперплоскость.

Мы получаем следующую интрпретацию метода наименьших квадратов:

  1. Первый геометрический факт Прогнозы, которые мы получаем по методу наименьших квадратов — это проекция исходного вектора зависимых переменных y на множество векторов, получаемых с помощью сложения с разными весами вектора из единичек, вектора x и вектора z.

  2. Второй геометрический факт Если я спроецировать вектор y на прямую, порождаемую вектором из единичек, и спроецировать y с крышечкой на ту же самую прямую, то эти проекции попадут в одну и ту же точку.

  3. \(TSS = ESS + RSS\)

  4. \(\frac{ESS}{TSS}= \frac{BC^2}{AB^2}=(\frac{BC}{AB})^2= (cos \varphi)^2\)

Геометрчиеские факты

Геометрчиеские факты

Коэффициент детерминации

Если в регрессию включён свободный член \(\beta_1\) то действуют следующие правила:

Правила при включении свободного члена в регрессию

Правила при включении свободного члена в регрессию

Эти правила позволяют придумать простой показатель качества работы модели — коэффициент детерминации.

Чем прогнозы точнее похожи на настоящие значения, тем меньше будут ошибки прогнозов и тем меньше будет сумма квадратов ошибок прогнозов RSS. Соответственно \(\frac{ESS}{TSS} = R^2\) будет примерно равен 1 если RSS будет у ноля. Соответственно мы получим коэффициент детерминации, который всегда лежит от 0 до 1.

Другими словами — коэффициент детерминации это доля объяснённого разброса в общем разбросе.

Правила при включении свободного члена в регрессию.

С одной стороны, коэффициент детерминации — это доля объясненной дисперсии игрек, доля объясненного разброса.

С другой стороны, коэффициент детерминации — это выборочная корреляция между прогнозами и настоящим игрек, взятое в квадрат.

Чем коэффициент детерминации выше, тем больше предсказание похожи на реальные значения. Чем коэффициент детерминации выше, тем выше доля объясненной дисперсии.

Пример с фертильностью. МНК в R

Посмотрим как зависит фертильность от других

t <- swiss
# нарисовать много диаграм рассеивания
# install.packages("GGally") # матрица диаграмм рассеяния
library("GGally")
## Loading required package: ggplot2
str(t)
## 'data.frame':    47 obs. of  6 variables:
##  $ Fertility       : num  80.2 83.1 92.5 85.8 76.9 76.1 83.8 92.4 82.4 82.9 ...
##  $ Agriculture     : num  17 45.1 39.7 36.5 43.5 35.3 70.2 67.8 53.3 45.2 ...
##  $ Examination     : int  15 6 5 12 17 9 16 14 12 16 ...
##  $ Education       : int  12 9 5 7 15 7 7 8 7 13 ...
##  $ Catholic        : num  9.96 84.84 93.4 33.77 5.16 ...
##  $ Infant.Mortality: num  22.2 22.2 20.2 20.3 20.6 26.6 23.6 24.9 21 24.4 ...
ggpairs(t)

Уже тут можно посмотреть много всяких корреляций.

Перейдём к оценке

model2 <-  lm(data = t, Fertility ~ Agriculture + Education + Catholic)
# КОэффициенты
coef(model2)
## (Intercept) Agriculture   Education    Catholic 
##  86.2250198  -0.2030377  -1.0721468   0.1452013
# Спрогнозированные значения
fitted(model2)
##   Courtelary     Delemont Franches-Mnt      Moutier   Neuveville 
##     71.35382     79.73758     86.36550     76.21257     62.05992 
##   Porrentruy        Broye        Glane      Gruyere       Sarine 
##     84.70365     77.94869     77.98965     82.07990     76.37831 
##      Veveyse        Aigle      Aubonne     Avenches     Cossonay 
##     81.01451     62.00804     65.34456     61.67811     67.20324 
##    Echallens     Grandson     Lausanne    La Vallee       Lavaux 
##     72.85406     71.22373     54.02437     62.00809     62.16632 
##       Morges       Moudon        Nyone         Orbe         Oron 
##     64.12130     72.47751     65.22299     69.41765     71.04507 
##      Payerne Paysd'enhaut        Rolle        Vevey      Yverdon 
##     66.61076     70.48740     64.27982     63.09324     68.48321 
##      Conthey    Entremont       Herens     Martigwy      Monthey 
##     81.11782     77.02791     80.38838     78.28372     84.09311 
##   St Maurice       Sierre         Sion       Boudry La Chauxdfnd 
##     75.54878     80.27332     73.53528     66.37864     74.87034 
##     Le Locle    Neuchatel   Val de Ruz ValdeTravers V. De Geneve 
##     70.52554     50.79966     71.80743     76.17918     35.30542 
##  Rive Droite  Rive Gauche 
##     52.99371     57.97821
# Остатки
residuals(model2)
##   Courtelary     Delemont Franches-Mnt      Moutier   Neuveville 
##    8.8461774    3.3624191    6.1345049    9.5874340   14.8400828 
##   Porrentruy        Broye        Glane      Gruyere       Sarine 
##   -8.6036475    5.8513084   14.4103471    0.3201012    6.5216935 
##      Veveyse        Aigle      Aubonne     Avenches     Cossonay 
##    6.0854871    2.0919628    1.5554442    7.2218873   -5.5032423 
##    Echallens     Grandson     Lausanne    La Vallee       Lavaux 
##   -4.5540632    0.4762715    1.6756342   -7.7080933    2.9336804 
##       Morges       Moudon        Nyone         Orbe         Oron 
##    1.3786987   -7.4775133   -8.6229883  -12.0176461    1.4549265 
##      Payerne Paysd'enhaut        Rolle        Vevey      Yverdon 
##    7.5892409    1.5125978   -3.7798150   -4.7932370   -3.0832083 
##      Conthey    Entremont       Herens     Martigwy      Monthey 
##   -5.6178154   -7.7279097   -3.0883806   -7.7837172   -4.6931098 
##   St Maurice       Sierre         Sion       Boudry La Chauxdfnd 
##  -10.5487835   11.9266828    5.7647206    4.0213575   -9.1703410 
##     Le Locle    Neuchatel   Val de Ruz ValdeTravers V. De Geneve 
##    2.1744592   13.6003353    5.7925740   -8.5791790   -0.3054172 
##  Rive Droite  Rive Gauche 
##   -8.2937095  -15.1782122
# ПОказатель RSS
deviance(model2)
## [1] 2567.884
# Показатель R квадрат
report <- summary(model2)
report$r.squared
## [1] 0.6422541
# Это тоже самое что
cor(t$Fertility, fitted(model2))^2
## [1] 0.6422541

Тесты

data("mtcars")

m1 <- lm(data = mtcars, mpg ~ disp + hp + wt)
summary(m1)$r.squared
## [1] 0.8268361
m2 <- lm(data = mtcars, mpg ~ cyl + hp + wt)
summary(m2)$r.squared
## [1] 0.84315
m3 <- lm(data = mtcars, mpg ~ disp + cyl + wt)
summary(m3)$r.squared
## [1] 0.832607
m4 <- lm(data = mtcars, mpg ~ disp + hp + cyl)
summary(m4)$r.squared
## [1] 0.7678877
m5 <- lm(data = mtcars, mpg ~ disp + hp + wt + am)
coef(m5)
##  (Intercept)         disp           hp           wt           am 
## 34.209443370  0.002489354 -0.039323213 -3.046747000  2.159270737
m6 <- lm(data = mtcars, mpg ~ disp + hp + wt)
summary(m6)$r.squared
## [1] 0.8268361
data("sleep")
var(tail(sleep, 11)$extra)
## [1] 3.618
max(sleep$extra) - min(sleep$extra)
## [1] 7.1
mean(sleep$extra)^3
## [1] 3.652264
sleep[8,]
##   extra group ID
## 8   0.8     1  8